Phase diagram of model anisotropic particles with octahedral symmetry 
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We computed the phase diagram for a system of model anisotropic particles with six attractive 
patches in an octahedral arrangement. We chose to study this model for a relatively narrow value 
of the patch width where the lowest-energy configuration of the system is a simple cubic crystal. At 
this value of the patch width, there is no stable vapour-liquid phase separation, and there are three 
other crystalline phases in addition to the simple cubic crystal that is most stable at low pressure. 
Firstly, at moderate pressures, it is more favourable to form a body-centred cubic crystal, which 
can be viewed as two interpenetrating, and almost non-interacting, simple cubic lattices. Secondly, 
at high pressures and low temperatures, an orientationally ordered face-centred cubic structure 
becomes favourable. Finally, at high temperatures a face-centred cubic plastic crystal is the most 
stable solid phase. 



I. INTRODUCTION 

Anisotropic patchy particles consisting of a repulsive 
core with some attractive sites have been used to study 
a variety of problems. The first anistropic patchy po- 
tentials were introduced as models of associative liquids 
(see, for example, Refs. HUHl). More recently, 
anisotropic patchy models have received renewed inter- 
est in the context of protein crystallization ,L£i£ i lo^wa 
Proteins are usually hard to crystallize, but, in order to 
determine the structure of a given protein, and hence 
its functionality, large crystals that can be used in high- 
resolution X-ray diffraction studies are needed.— In this 
respect, theoretical studies aiming to predict the condi- 
tions which favour crystallization would be very valuable. 
Even though some significant progress has already been 
made using simple isotropic potentials , 15 ! 16 it is known 
that protein interactions are short-ranged and highly 
anisotropic. For example, most protein crystals have 
packing fractions that are much lower than the close- 
packed solids typically favoured by isotropic potentials^ 
So far, studies using anisotropic models have shown 
that the fluid-fluid coexistence moves to lower temper- 
atures as the interactions become more anisotropic (ei- 
ther by decreasing the number of patches or making them 
smaller) 1 ^ and can even become metastable.— The intro- 
duction of anisotropy can also induce the stability of mul- 
tiple solid phases, including orientationally ordered and 
plastic phases* 1 - 1 - 

Anisotropic interactions have also received much at- 
tention in the context of the development of new ma- 
terials. In particular, there has been increasing inter- 
est in the fabrication of colloids 1 ^^ 2 ^ 2 ^ 2 ^^ and 
nanoparticle o 27 ! 28 ' 29 with anisotropic interactions, one of 
the goals being to tailor their interactions so that they 



are able to self-assemble into a given target structure. 
One set of targets that has received particular attention 
is low-density colloidal crystals, e.g. colloidal diamond, 
because of their potential as photonic materials^ This 
work has stimulated a number of recent theoretical stud- 
ies that have begun to address the question of how such 
patchy interactions can be used to control the crystalliza- 
tion behaviou r 31 ! 32 and the self-assembly of finite objects 
of given size and symmetryj 33 i 34 i 35 i 36 Some of the latter 
studies have also been motivated by a desire to under- 
stand biological self assembly, such as the formation of 
virus capsids. 

The final area where anisotropic patchy models have 
been the subject of recent interest is in the study of the 
dynamics of supercooled liquids j 37 ' 38 i 39 ' 40 In particular, 
as patchy models tend to move the liquid-vapour coex- 
istence line to lower temperature, and to lower packing 
fractions,— they make it much easier to study gels, i.e. 
dynamically-arrested states that have low density and 
where the arrest is due to formation of energetically sta- 
ble bonds, rather than caging of the atoms in a densely- 
packed environment. 

In spite of these studies, there is still much to be learnt 
from a fundamental point of view about how anisotropic 
interactions affect the thermodynamics and dynamics of 
a system. In particular, it is not fully understood how 
the arrangement of the patches will affect the phase dia- 
gram. It seems reasonable that if the patches are located 
at positions that favour the local environment of a given 
crystalline structure, crystallization into that structure 
would be favoured. On the other hand, a random dis- 
tribution of the patches, as perhaps may be the case for 
the surface of a protein, will normally lead to a situa- 
tion where the local order is not compatible with any 
crystalline lattice and, therefore, crystallization will be 
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hindered. However, relatively little is known about how 
the specificity of the angular interactions will affect the 
relative stability of the phases or their accessibility due 
to kinetic effects. 

In this paper we study a model of patchy particles that 
we have previously used to study the self-assembly of 
monodisperse clusters, 3 ? and the kinetics of crystalliza- 
tion in two and three dimensions^ One of the intriguing 
results of the latter work was that crystallization of parti- 
cles with six octahedrally-arranged patches into a simple 
cubic (sc) crystal appeared to be much easier than crys- 
tallization of particles with four tetrahedrally-arranged 
patches into a diamond lattice. Our hypothesis is that 
this difference in crystallization kinetics reflects the ab- 
sence of frustration in the octahedral system, whereas 
in the tetrahedral system the preferred local order differs 
from the global crystalline order, thus frustrating crystal- 
lization. In order to understand this further a systematic 
study of the nucleation behaviour for these two systems is 
required, and a necessary precursor for such work is the 
computation of the phase diagram. This is one of the 
motivations for the present study, where we compute the 
phase diagram for particles with an octahedral arrange- 
ment of the patches. Furthermore, the phase diagram 
will also be of considerable interest in its own right, with 
the potential for competition between lower- and higher- 
density crystalline phases, and the possibility of plastic 
phases that have translational, but no orientational or- 
der. 



II. METHOD 
A. Model 

Our model consists of spherical particles with a given 
number of attractive patches whose geometry is speci- 
fied by a set of patch vectors. The total potential can 
be written as a sum of two-body terms that depend on 
the distance between two particles , but also on their 



relative orientations Jl; and fi,-: 
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The interaction between two particles is described by a 
potential with an isotropic repulsive core and an angular- 
dependent attractive term: 



V(ry-,f2i,%) = 
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where Vlj(t) is the Lennard- Jones potential 
V LJ (r) =4e 
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e is the pair well depth and Vlj(<Jlj) — 0. Additionally, 
for computational efficiency, the potential is truncated 




FIG. 1: A schematic representation of the geometry of the 
interaction between two particles. For clarity, we depict a two- 
dimensional analogue of the three-dimensional model used in 
this work. In this two-dimensional model, the particles have 
four patches arranged regularly with their directions described 
by the patch vectors, pi. In the particular case shown in the 
figure, patch 4 on particle i interacts with patch 2 in particle 
j because they are the closest to the interparticle vector. 



and shifted using a cutoff distance of 2.5 <jlj- The attrac- 
tive interaction is modulated by a product of Gaussian 
functions that are centred at the position of each patch: 



fli, flj) = exp 



2a 2 



exp 



2ct 2 
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where a is the standard deviation of the Gaussian, 9k,ij 
(Oi.ji) is the angle formed between patch k (I) on atom i 
(j) and the interparticle vector ( Yji), and k m i n (l m in) 
is the patch that minimizes the magnitude of this angle. 
The interaction is a maximum when both patches are 
pointing at each other along the interparticle vector 
and it will decrease as the particles deviate further from 
this equilibrium orientation. A schematic representation 
of two such interacting particles is provided in Fig. [1] 

The angular dependence of Eq. |4] mimics the orienta- 
tional dependence that exists in short range directional 
forces as is the case for hydrogen bonding. The use of 
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means that for a given pair of par- 



ticles only a single patch on each particle is involved in 
the interaction, i.e. the possibility of 'double hydrogen 
bonding' is removed by the use of Eq. |4j 

One of the advantages of this potential is its simplic- 
ity. It is easy to implement and it is computationally 
not very expensive to evaluate. Furthermore, the model 
includes the anisotropy in a very simple and flexible way. 
Simply by changing the number and position of patches, 
each of which is defined by a vector in the particle ref- 
erence system, it is possible to obtain a wide range of 
anisotropic potentials. Another advantage is that the 
well-characterized Lennard-Jones potential can be ob- 
tained as a limiting case of the current model when the 
width of the patches becomes increasingly large. This 
feature can be particularly useful if one wants to study 



the effect of the anisotropy, going from a very anisotropic 
model to the isotropic limit. 

We shall use reduced units throughout, so that U* — 
U/e, p* = p/(e/alj), T* = T/(e/k B ) and p* = pa 3 LJ . 
Consequently, the only parameter that needs to be spec- 
ified to fully characterize the interaction between two par- 
ticles is a, i.e., the angular width of the attractive patches 
(see Eq. 0J. In the present calculations, we use particles 
with six patches that have an octahedral arrangement, 
and since we wish to study the model in a regime where 
there is a stable low-density crystal, we have chosen to 
use a relatively narrow patch width, namely a =0.3 radi- 
ans. In previous work^ it has been shown that, for this 
value of a and at not very high pressures, a low density 
sc crystal is formed spontaneously from the liquid at suf- 
ficiently low temperature. In particular, at a constant 
pressure p* =0.1 crystallization occurred at T* «0.17. 



B. Solid structures 

There are several crystalline structures that one might 
expect to be stable for particles with an octahedral ar- 
rangement of the patches. The most simple structure is 
a simple cubic crystal. In this structure, each of the six 
patches points directly at one of the six nearest neigh- 
bours, and none of the interactions will be frustrated. 
However, this crystal has a relatively low density, e.g. 
if the nearest-neighbour separation is equal to the min- 
imum in the pair potential, p* — l/v2 = 0.707, and 
so it is expected that new denser crystalline phases will 
appear when the system is exposed to moderate to high 
pressures. 

Two possible higher-density crystals are the body- 
centred-cubic (bcc) and face-centred-cubic (fee) lattices. 
However, in these structures, due to the symmetry of 
the potential, it is not possible to orient the particles so 
that each patch is pointing directly towards one of its 
nearest neighbours. For a bcc crystal, it is possible to si- 
multaneously align only two of a particle's patches with 
its first neighbours. However, an orientationally ordered 
structure can be formed when the six patches of a given 
particle are pointing towards its six second nearest neigh- 
bours (see Figure [2]). This structure can also be viewed 
as two interpenetrating simple cubic lattices that, as we 
will show, almost do not interact with each other. As 
the patches form an angle of ir/4 radians with the first 
neighbours, there will almost be no interactions between 
first neighbours, except when they are closer than the 
repulsive barrier (ctlj). Furthermore, when the nearest- 
neighbours separation is <tlj, the distance betweeen sec- 
ond nearest neighbours is 2oLjj\fZ — \.15A7aLj which 
is only slightly longer than the minimum in the pair po- 
tential that occurs at 2 1 / 6 <7lj = 1.1225ctlj. Therefore, 
the energetic cost of interpenetrating the two lattices is 
relatively small, but as the density is significantly higher 
(p* = 1.299 for the above geometry) the bcc lattice will 
have a more favourable enthalpy at moderate pressures. 
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FIG. 2: (Colour online) Unit cells of the (a) bcc and (b) 
orientationally ordered fec-o structures. As mentioned in the 
text, in both cases the patches are aligned with the second 
neighbours. 



In some senses, this structure has similarities to the struc- 
tures of ice VII and ice VIII. Both consist of two inter- 
penetrating cubic ice lattices that do no have any inter- 
connecting hydrogen bonds between them4I 

For the fee lattice, it is also possible to generate an 
orientationally-ordered structure, where the six patches 
are pointing towards six of the second nearest neighbours 
(we refer to this structure, which is depicted in Figure 
as fec-o). However, it is significantly higher in energy 
than the bcc lattice, because in this case the ratio of the 
second to the first nearest-neighbour distance is \/2, and 
so if there is to be no repulsion between nearest neigh- 
bours (e.g. if their separation is a^j and p* = then 
the second nearest neighbours are significantly further 
apart than the minimum in the pair potential. 

It is also possible to generate an orientationally- 
ordered tetragonal crystal, in which the centres of mass 
of the particles are disposed as in a fee lattice, but where 
four of the patches point towards first neighbours and 
the other two patches point towards second neighbours. 
This structure has a body-centred tetragonal unit cell 
(a = b =/= c, where a, b and c are the moduli of the lat- 
tice vectors) , as the cubic symmetry of the fee lattice has 
been broken due to the orientation of the particles (see 
Figure 1 in Reference for an explanation of how to get 
the tetragonal unit cell in an analogous structure for a 
model of oppositely charged colloids). However, in NpT 
simulations in which each of the edges of the box was 
allowed to vary independently, we found the tetragonal 
crystal to undergo a transformation to an orientationally 
ordered bcc lattice. This type of deformation has been 
previously observed in other systems, and it is an ex- 
ample of a martensitic transition.™ The transformation 
occurs by a shortening of the c edge, until c = a and 
cubic symmetry is recovered. In view of these results, we 
have not further considered the tetragonal structure in 
our calculations. Nevertheless, one cannot exclude the 
possibility that there might be some range of pressure 
and temperature where the tetragonal structure is ther- 
modynamically stable. 
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At sufficiently high temperatures when the attractive 
interactions become negligible, a translational ordered, 
but orientationally disordered fee structure, i.e. a plastic 
phase, will also appear. We refer to this structure as 
fec-d (or PC). 



C. Equation of state for the fluid and solid phases 

The equation of state for the solid and the liquid phases 
was calculated using Monte Carlo (MC) simulations in 
the NpT ensemble. Between 10 and 20 simulations were 
performed to determine the equation of state along each 
isotherm. Each simulation consisted of 40000 MC cycles, 
following an equilibration period of the same length. A 
MC cycle was defined as N attempts to translate a par- 
ticle, plus N attempts to rotate a particle and two at- 
tempts to change the volume of the simulation box, N 
being the number of particles in the system. During the 
equilibration period, the maximum translational and ro- 
tational displacements were adjusted to obtain a 40% 
acceptance probability and the maximum volume change 
was chosen to obtain a 30% acceptance probability. In all 
the simulations, the output configuration of a given state 
was used as the input for the following simulation. The 
number of particles used in our calculations was 216 for 
the sc phase, 250 for the bec phase, 256 for the fec-o and 
fec-d phases and 250 for the fluid phase. In all the cases, 
the length of the box was larger than twice the cutoff in 
the potential. Since all the solid structures that we con- 
sidered were cubic, we performed NpT simulations with 
isotropic scaling. 



D. Free energy calculations 

At at given temperature, coexistence between two 
phases occurs at the pressure p where the chemical po- 
tential p is the same for both phases. Since p/k B T = 
G/Nk B T = A/Nk B T + pV/Nk B T and the compressibil- 
ity z = pV / Nk B T can be obtained via NpT simulations, 
a procedure to determine the Helmholtz free energy is 
needed. 

The free energy of the fluid phase was estimated by 
integration from a very low density state, where the fluid 
can be considered to behave as an ideal gas: 



A(p) _A«{p) 



Nk B T Nk B T 



1 



dp 



(5) 



where A id {p)/Nk B T = bx(pA 3 ) - 1. We took A = a L j 
because its value does not affect the coexistence proper- 
ties. 

For the solid phases, we used the method described 
by Frenkel and Lad d 43 i 44 , as extended to non-spherical 
potentialsi^ ' 45 i 46 ' 47 i 48 In this method, the free energy is 
obtained by integration from the interacting Einstein 



crystal, whose interactions are described by the Hamil- 
tonian: 



H(X{,X* 2 ) H 
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On the right hand side of this equation, H is the orig- 
inal potential and the second and third terms are, re- 
spectively, an orientational and a translational field that 
tend to keep the particles at the positions and orienta- 
tions of the perfect orientationally ordered lattice. In the 
term for the orientational field, \l/ a is the minimum an- 
gle formed by any of the vectors that define the position 
of the patches in the particle's reference system with re- 
spect to the x axis of a fixed reference system and ^f, is 
the analogous quantity with respect to the y axis, where 
the fixed reference system has been chosen to be coin- 
cident with the orientation of the patches in the perfect 
lattice. As the patch that defines both angles must not 
be the same, when the same patch yields the minimum 
angle with both the x and y axes, the patches that form 
the second least angles were also computed. The patches 
that will contribute to the orientational field will be the 
pair that leads to the lowest energy for this term. Note 
that the orientational field has been chosen so that it has 
the same symmetry as the particles. 

In the term for the translational field, ro,j is the lattice 
position of particle i (as given by its centre of mass). 
A^ = \\jk B T and A?; = X 2 /(k B Ta LJ ) are the cou- 
pling parameters of the translational and orientational 
field, respectively. For practical reasons, we chose to use 
X'l = X'2 = A*, i.e. both fields are switched on simulta- 
neously. With this choice, the integral to the reference 
translational and orientational Einstein crystals can be 
performed at the same time, thus reducing the number 
of simulations needed to perform the numerical integra- 
tion to the reference system. However, other choices of 
Ai and A2 are equally valid. Therefore, by using the same 
coupling parameter for both fields, Eq. [S]can be simpli- 
fied into the following form: 
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The free energy difference between our model and the 
reference Einstein crystal (AA2) is given by: 



&A 2 r° / d(H(X*)/Nk B T) 
Nk B T J x . \ dX* 



dX* 



(8) 



This integral was evaluated numerically using a Gauss- 
Legendre quadrature formula, with 10 or 20 points de- 
pending on the case (a larger number of points is needed 
for the orientationally disordered plastic phase). 
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The free energy of a non-interacting Einstein crystal is 
already known^ and the free energy difference between 
an interacting and a non- interacting crystal (AAi) can 
be calculated numerically by averaging the interparticle 
energy over a simulation of the non-interacting Einstein 
crystal^&is 



Nk B T 



exp 



U 



Nk B T 

JV-l JV 

E E 
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k B T 
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where Uq is energy of the perfect lattice. The maximum 
value of A* in Eq. [5] was chosen so that the structure 
of the interacting Einstein crystal defined by Eq. [5] was 
very similar to the perfect orientationally ordered lattice. 

The free energy of the orientational held can be esti- 
mated numerically by integrating its partition function 
over all the orientations^^ 



i 5 




FIG. 3: Determination of the coexistence point between the 
fluid (solid line) and the bcc (dashed line) phases at T* = 
0.243. 



Nk B T 



In j exp {-A* [sin(* ) 2 + sin(* 6 ) 2 ] } x 

sm(a)dad(f>dj) (10) 



where \l/ a and 'J'b depend on the three Euler angles, a, 
(j) and 7. This integral was evaluated numerically using 
the Monte Carlo integration method and using at least 
10 9 points. 

The final expression for the free energy of the solid is: 



.4' 
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Einstein 



AAi + AA 2 + AA 3 



(11) 



where A Elnstein is the sum of the free energy of the trans- 
lational Einstein crystal— plus the free energy of the ori- 
entational held A onent . The term A^3 accounts for the 
fact that the integration to the Einstein crystal was eval- 
uated by performing simulations in which the centre of 
mass of the system was fixed (see Ref. |43T ). 

Once the free energy is known for a given state, the free 
energy along an isotherm can be obtained by integrating 
the equation of state: 



A(p 2 ) A( Pl ) 
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(12) 



The coexistence points are determined by imposing the 
conditions of equal pressure and chemical potential, 
and can be obtained by plotting the chemical poten- 
tial against the pressure for both phases. Coexistence 
is where the curves for the two phases cross, and Fig. 
[3] shows such a plot for the bcc and the fluid phases at 
T* = 0.243. 



E. Coexistence lines 

The coexistence lines have been located using 
the Gibbs-Duhem integration method introduced by 



Kofke i 50 ! 51 In this method the coexistence line is obtained 
by integration of the Clausius-Clapeyron equation: 



dp \ / Au + pAv 
dTJ = { T~Av 



(13) 



where Am and Av are the molar energy change and molar 
volume change, respectively, between the two coexisting 
phases. 

The Clausius-Clapeyron differential equation was 
solved numerically with a fourth order Runge-Kutta al- 
gorithm, and using as the initial conditions the coexis- 
tence points obtained with the thermodynamic integra- 
tion method. The value of the integrand in each of the 
four steps of the Runge-Kutta method was evaluated by 
means of NpT simulations that consisted of 10000 MC 
cycles, following 10000 cycles of equilibration. 



F. Direct coexistence simulations 

As a check of the above calculations, the melting points 
of the three solid phases were also estimated using the 
direct coexistence method, hrst proposed by Ladd and 
Woodcock.™ In this work, we will follow the same pro- 
cedure as the one described in Rcf. 53. A crystal with 
around 400-500 particles was generated and subsequently 
equilibrated in the iVpT-ensemble at a given pressure and 
temperature. This configuration was copied and heated 
to obtain a liquid configuration, which was then equi- 
librated at the same conditions as the solid but using 
a A^T-ensemble in which the volume of the simulation 
box can only change by modifying one of the box lengths, 
which we chose to be the box length along the z direc- 
tion. Therefore, both the solid and liquid phases have 
the same periodic conditions along the x and y axes, and 
a solid-liquid interface can be built by simply joining the 
liquid and solid configurations along the x — y plane. 
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FIG. 4: (Colour online) The variation of the internal energy 
per particle (w* = U* /N) during NpT simulations of a box 
containing the bcc solid in contact with the fluid at p* = 1.0. 
A few trajectories at different temperatures are shown. 




is below the melting, the solid phase will act as a nucle- 
ation seed and the crystal will grow at the expense of the 
liquid, resulting in a lowering of the system's energy. If 
it happens that the temperature is equal to the melting 
temperature, then the solid would be in equilibrium with 
the liquid, and the energy would remain constant. Using 
this procedure, we can determine a temperature inter- 
val which give upper and lower bounds for the melting 
temperature. 

The results of NpT simulations of a box containing 
the bcc solid and liquid phases at p* — 1.0 are shown in 
Figure 2J In this Figure, it can be seen that the melt- 
ing temperature is approximately T^ elt = 0.243. At this 
temperature the internal energy oscillates around an av- 
erage value. At temperatures higher than T* = 0.243, 
the energy increases until it reaches a constant value, 
that corresponds to the situation when all the solid has 
melted. Finally, at temperatures lower than T* = 0.243, 
the energy decreases, again until it reaches a constant 
value, which, in this case, means that all the fluid has 
frozen. It is worth noting that both the time for the 
fluid to freeze and the time for the solid to melt become 
shorter as the temperature of the system moves further 
away from the melting temperature. Figure [5] shows the 
initial configuration of the simulation box, as well as two 
final states, one above and one below the melting. A 
visual inspection shows that the final configuration at 
T* =0.245 corresponds to an homogeneous fluid phase 
(Fig. [5] (b)) and, on the contrary, at T* =0.240 it corre- 
sponds to a perfect crystal (Fig. [5] (c)). 



III. RESULTS 




FIG. 5: (Colour online) (a) Snapshot of the initial configu- 
ration of the simulation box containing the bcc and the fluid 
phases in contact, (b) and (c) Snapshots of the final con- 
figurations for T* =0.245 and T* =0.240, respectively. The 
pressure was set to p* =0.1. 



For a given pressure, the melting temperature can be 
estimated performing NpT MC simulations at different 
temperatures. If the temperature is above the melting 
point, the solid phase will melt and the energy of the 
system will increase. On the contrary, if the temperature 



Let us start by presenting the results for the fluid 
phase. The free energy of the fluid phase was obtained 
by thermodynamic integration from the very low density 
limit, where the fluid can be considered to behave as an 
ideal gas, to high densities (Equation [5]) . The integrand 
of this equation was evaluated by NpT simulations for 
different densities, and the results were fitted to a poly- 
nomial of degree six: 



P* 



= a + aip 



*2 , *3 i *4 

0-2P + a 3 p + Q4P 



(14) 



a 5 p + a e p 



The coefficients resulting from this fit for different 
isotherms are shown in Table |TJ We have found no evi- 
dence of a vapour-liquid transition in any of the isotherms 
studied, even at the lowest temperature we studied, T* — 
0.200. The coefficient ao provides an estimation of the 



second virial coefficient of the model (£?2 = Urn 



z-l\ 

p^o —)■ 

As it can be seen, the Boyle temperature, i.e., the tem- 
perature for which the second virial coefficient vanishes, 
is located between 0.200 and 0.243. A more precise es- 
timate can be obtained by determining B2 from numeri- 
cal integration. Using this procedure we have computed 
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T* 


do 


ffli 


a 2 


Ct3 


0,4 


a 5 


0,6 


0.200 


-0.6133262 


2.3616349 


9.5903228 


-55.187794 


138.527889 


-122.51950 


-399.33624 


0.243 


0.1538809 


2.6642095 


-2.2559535 


2.8910917 


29.6712406 


-39.785495 


18.816227 


0.500 


1.2862246 


1.6415799 


8.4391218 


-30.150196 


74.267261 


72.615067 


28.555569 


3.00 


1.6574717 


2.3183855 


-1.5757928 


10.982981 


-14.129931 


10.595277 


-2.4231455 



TABLE I: Coefficients obtained by fitting the integrand of Equation ([5| to a polynomial of degree six (Eq. (JT4J). 



1.00 



0.75 



*h-0.50 



0.25 



O.Ofl 



o a=0.7rad. 
* a=0.5 rad. 



— o 
o 

o 



* 

-* 



po 



°o, 



00 



0.25 



0.50 
P 



0.75 



1.00 



FIG. 6: Dependence of the vapour-liquid coexistence curve 
on the patch width a. At lower a we were unable to find the 
(metastable) coexistence curves. 



B 2 at T* = 0.200, 0.220, 0.230, 0.240, 0.500, obtain- 
ing B 2 /<tIj = -0.556,-0.132,0.031,0.172,1.386. From 
these results, we estimated that T^ oyle w 0.228. An ap- 
proximation to the critical temperature (T c ) can be ob- 
tained from the Boyle temperature. It is known that for 
another associative model, the primitive model of water, 
TBoyie/T c ~ 1.244^°. Even though it is not clear how 
the anisotropy of the model will affect this quotient, we 
can obtain a rough estimate of the critical temperature. 
Taking T* oyle w 0.228, T* « 0.18. This estimate is con- 
sistent with the fact that the equation of state of the fluid 
phase along the isotherm T* = 0.200 docs not show any 
indication of a vapour-liquid equilibrium. 

The possibility of a vapour-liquid phase separation 
was further explored by means of Gibbs ensemble (GE) 
simulations! 54 ' 55 ' 56 Our GE simulations consisted on 5 x 
10 5 MC cycles for equilibration plus 10 6 MC cycles for 
obtaining averages, where a MC cycle was defined as 
N attempts to translate a particle, N attempts to ro- 
tate a particle, two attempts to change the volume and 
about 200 attempts to exchange particles between the 
two boxes. We simulated a system with 256 particles at 
a constant total number density p* = 0.3. Several patch 
widths (<r =0.3, 0.5 and 0.7 radians) were considered and 
we found that, consistent with previous calculations j 7 -* 1 ^ 
the vapour-liquid coexistence curve is very sensitive to 
the width of the patches, and that the critical temper- 



ature, which was estimated by fitting the GE results to 
the law of rectilinear diameter and to the critical expo- 
nent scaling law ; 57 ' 58 decreases as the patches become 
narrower (see Fig. [5]). For the particular case a = 0.3 
radians, we were unable to locate any vapour-liquid co- 
existence points, because, in this case, the simulations 
need to be performed at very low temperatures, where it 
is particularly difficult to obtain well-converged results. 
In any case, the observed trends indicate that the criti- 
cal temperature for a = 0.3 must be very low and will 
fall below the fluid-sc coexistence curve, which is greater 
than T* = 0.2 up to fairly low densities. These results 
suggest that if vapour-liquid phase separation does exist, 
it must be metastable with respect to freezing. This is 
a relevant finding, because several experimental studies 
have shown that globular proteins also exhibit a vapour- 
liquid phase separation that is metastable with respect 
to solidification ! 59 ' 60 Moreover, ten Wolde and Frenkel 
have suggested that crystallization occurs more rapidly 
in the proximities of this metastable critical point J£ In 
that sense, it would be interesting to determine what is 
the maximum patch width for which the phase separa- 
tion between a low density fluid and a high density fluid 
is metastable. 

The free energies of all the solid phases at several ther- 
modynamic states are given in Table [TTJ In some cases, 
we computed the free energy for two states on the same 
isotherm, as these can then be used to test the thermo- 
dynamic consistency of our results (i.e. the free energy 
difference between two states as calculated using the Ein- 



Structure 


T* 


P* 


^max 


Atot 


AA 2 




U 


Sc- 


0.100 


0.762 


25000 


-13.417 


-12.644 


-29.269 


-29.352 


heie 


0.100 


1.230 


20000 


-11.637 


-10.784 


-28.667 


-28.723 


fec-o 


0.100 


1.360 


20000 


4.355 


-11.689 


-12.455 


-12.508 


se 


0.200 


0.763 


20000 


-1.095 


-14.284 


-14.616 


-14.669 


bec 


0.200 


1.175 


20000 


0.340 


-13.754 


-13.719 


-13.740 


bec 


0.200 


1.210 


20000 


0.690 


-13.001 


-14.122 


-14.146 


bec 


0.243 


1.147 


20000 


2.055 


-14.754 


-11.004 


-11.020 


fec-d (PC) 0.500 


1.204 


20000 


5.452 


-20.326 


-1.971 


-1.971 


fec-d (PC) 


3.00 


1.283 


20000 


5.488 


-21.890 


-0.372 


-0.373 


fec-d (PC) 


3.00 


1.376 


20000 


6.515 


-20.941 


-0.295 


-0.313 



TABLE II: Free energy of the solid phases, as obtained by 
thermodynamic integration from the Einstein crystal. The 
free energies (Atot, AA 2 and AAi) and the lattice energy 
(Z7o) are given in units of NksT. 
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Solid phase T* bo bi 



sc 0.100 26.7846382 -69.2645290 6.67906462 54.0334053 

bcc 0.100 612391.025 -2010382.65 2475295.23 -1354791.73 278126.549 

fcc-o 0.100 -631690.341 1383537.04 -1010190.80 245893.687 

bcc 0.243 -570.084541 1657.42401 -1613.44221 525.947816 

fcc-d (PC) 3.00 -1410.40288 3664.08695 -3179.10036 943.459700 

TABLE III: Coefficients of the polynomial fit to the equation of state of the solid phases. The points were fitted to a third-degree 
polynomial, except for the bcc structure at T* = 0.1, for which a fourth-degree polynomial significantly improves the fit. 



stein crystal must be the same as that obtained from in- 
tegration of the equation of state). In particular, for the 
bcc phase at T* = 0.2, the difference in free energy be- 
tween the states at reduced densities 1.210 and 1.175 is 
0.3507VfcsT as calculated using the Einstein crystal ap- 
proach, which compares well with the value obtained by 
integrating the equation of state (0.355-/VfcsT). For the 
plastic fcc-d phase at T* = 3.0, the agreement is also 
good (the free energy difference between the states at re- 
duced densities 1.376 and 1.283 computed using the two 
methods is 1.027 Nk B T and 1.034 Nk B T). This ther- 
modynamic consistency provides positive confirmation of 
the reliability of the calculations. 

The free energy of the solid phases can be obtained 
at any other point on the isotherm using Eq. 1121 The 
equations of state (p* (p*)) for the solid phases were ob- 
tained by performing NpT simulations at different ther- 
modynamic states along a given isotherm and fitting the 
results of the simulations to a polynomial of the form: 

p*(p*) = 6 + b lP * + b 2 p* 2 + b 3 p* 3 + b 4P * 4 . (15) 

The resulting data was fitted to a polynomial of degree 
three, except for the bcc structure at T* = 0.1, for which 
a fourth-degree polynomial lead to a much better fit. 
The coefficients resulting from this fitting procedure are 
shown in Table HTT1 

Using these free energies, we calculated the coexistence 
points between all the phases that can be at equilibrium, 
and the results are shown in Table ITVl In some cases, we 



Phase 1 


Phase 2 


T* 


P* 


Pi 


P*2 




fluid 


fcc-d (PC) 


3.00 


48.0 


1.182 


1.270 


* 


fluid 


fcc-d (PC) 0.500 


5.55 


0.975 


1.070 




fluid 


bcc 


0.243 


1.013 


0.767 


1.107 


* 


fluid 


bcc 


0.200 


0.339 


0.615 


1.123 




fluid 


sc 


0.200 


0.048 


0.228 


0.672 




sc 


bcc 


0.200 


0.709 


0.711 


1.136 




sc 


bcc 


0.100 


0.444 


0.716 


1.208 


* 


bcc 


fcc-o 


0.100 


26.66 


1.119 


1.407 


* 



TABLE IV: Coexistence points obtained using the thermo- 
dynamic integration method. The points marked with an as- 
terisk were used as the starting points for the Gibbs-Duhem 
approach. The other points served to test our calculations. 



have calculated the coexistence between two phases at 
two different temperatures, in order to verify that Gibbs- 
Duhem integration was able to give accurate results even 
for regions quite far from the starting point. 

As a further test of our calculations, the melting points 
of the solids have also been calculated using the direct 
coexistence method. The melting points for all the solid 
phases obtained with this technique are shown in Table 
fVl The agreement between both methods is fairly good, 
the differences being of the order of 1%. 

Using the coexistence points in Table IIVI as initial 
conditions, we traced the coexistence curves with the 
Gibbs-Duhem method. Although we usually integrate 
the Clapeyron equation, as given by Eq. [131 we some- 
times found it more convenient to integrate the equation 
dT/dp = (TAv)/(Au + pAv). Some of the points ob- 
tained with the Gibbs-Duhem are shown in Table IVll 

The T — p and p — T phase diagrams are shown in 
Figures [7] and respectively. The dashed line in the 
diagrams shows a transition between the orientationally- 
ordered and disordered fee structures. These points have 
been estimated by heating the ordered fcc-o structure and 
monitoring the internal energy. This quantity exhibits a 
quite abrupt change when orientational order is lost, and 
the transition temperature was chosen as the tempera- 
ture where the internal energy curve shows an inflection 
point. 

As expected, the phase diagram shows multiple solid 
phases. Firstly, at high temperature, where the be- 
haviour is dominated by the repulsions between the par- 
ticles, the fluid freezes into a plastic crystal phase, the 
fcc-d phase (i.e., a lattice with fee structure with respect 
to the centre of mass, but with orientational disorder). 
Secondly, at intermediate temperature, it freezes into a 



Solid 


Direct coexistence 


Free energy calculations 




p* T* 


V* 


T* 


sc 


0.758 0.232 ± 0.004 


0.758 


0.229 


bcc 


1.000 0.243 ± 0.004 


1.013 


0.243 


fcc-d (PC) 


48.0 2.97 ± 0.02 


48.0 


3.00 



TABLE V: Melting points obtained using the direct coexis- 
tence method. For comparison, the coexistence points ob- 
tained from free energy calculations (see table lIVp are also 
shown. Note that the sc-fluid coexistence point was obtained 
using the Gibbs-Duhem method (see table [VlJ) . 
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L IlclSC L 


Phase 2 


p* T* 


* 
Pi 


* 


f>2 


u 2 


fluid 


fcc-d (PC) 48.00 2.500 


1 149 


3 Qfl1 


1 937 
1 .ZD I 


3 1 QR 

o. iyo 


fluid 


fcc-d (PC) 38.19 2.000 


1 112 


i . y lo 


1 907 
1 . ZU I 


1 . 440U 


nuid 


fcc-d (PC) 12.65 1.000 


1 031 


n 3^zi 

U.O04 


1 1 9f^ 
1 . 1ZO 


n 1 70 


fluid 


fcc-d (PC) 0.624 0.500 


-0.206 


-0.482 


1.057 


-0.303 


bcc 


fcc-d (PC) 


3.42 0.336 


1 19^ 


1 R39 


i oa^ 


O A70 




fcc-d (PC) 


5.00 0.355 


1 Kfi 


- l.UOU 


1 121 


41 7 


bcc 


fcc-d (PC) 


7.50 0.363 


1 1 CM 

±. iy4t 




1 1 Qf; 

i . iyo 


O 3^3 
-U.oOo 




fcc-d (PC) 10.25 0.365 


1 224 


-l.Olo 


1 249 


n 981 

-VJ.ZU1 


bcc 


fcc-d (PC) 15.58 0.325 


1 9f^ 
l.ZOO 


1 so^ 


1 399 

i .ozz 


n 1 1 q 
-u. i iy 


bcc 


fcc-d (PC) 22.36 0.250 


1 9CM 

i.zy^t 


1 CMQ 


1 374 
1 . I 44 


n oas 




fcc-o 


26.59 0.040 


1.303 


-1.648 


1.410 


-0.353 


nuid 


bcc 


2.915 0.320 


Q97 


n 3Q0 


1 191 
1 . 1Z1 


-1 .o*±\j 


nuid 


bcc 


2.074 0.290 


o 873 


n A99 


1 113 
1 . 1 lo 


1 7^Q 

-i . ( oy 


fluid 


bcc 


0.778 0.230 


n 79Q 

u. / zy 


o ZLQA 


1 110 


-i .yu^t 


sc 


bcc 


0.751 0.225 


n 7DQ 

u. / uy 


9 1 X.A 
-Z. 104 


1 1 1 Q 
1 . 1 lo 


9 nn/i 

-Z.UU4 


sc 


bcc 


0.582 0.150 


0.714 


-2.493 


1.173 


-2.381 


sc 


bcc 


0.382 0.080 


0.716 


-2.746 


1.222 


-2.677 


fluid 


sc 


0.758 0.229 


0.728 


-0.494 


0.709 


-2.133 


fluid 


sc 


0.200 0.218 


0.486 


-0.378 


0.677 


-2.165 


fluid 


sc 


0.010 0.181 


0.059 


-0.070 


0.675 


-2.333 



TABLE VI: Some of the coexistence points obtained with the 
Gibbs-Duhem method. 



bcc structure. Finally, at low temperature, it freezes into 
a low density sc solid. The sc structure is only stable 
at fairly low pressures, as it is possible to obtain a more 
dense phase, the bcc crystal, just by introducing a single 
atom at the centre of the unit cell without a large ener- 
getic penalty (see the discussion above and Table HI)) . The 
bcc phase remains stable up to considerably higher pres- 
sures. As the bcc structure is energetically much more 
favourable than the fcc-o structure (Table [TTJ) , the bcc 
crystal is only destabilized at densities for which the first 
neighbours are at distances close to the LJ repulsive core 
<jlj- In the fcc-o structure, the patches are pointing to 
the second neighbours, which are at a distance consider- 
ably larger than the LJ minimum (approximately v^ctlj 
or larger). 

It is worth noting that the sc, bcc and fcc-o phases are 
all stable at T* — 0. At this temperature, the sc and bcc 
structures are both stable over a finite, but very small 
range of density, because under these circumstances the 
solid becomes almost incompressible. 

The phase diagram exhibits at least two triple points 
whose thermodynamic states are given in table IVIII (we 
have not studied in detail the triple point where the bcc, 
fcc-d and fcc-o phases coexist). At one of these triple 
points, the sc crystal coexists both with the liquid and 
the bcc solid phases. This triple point is somewhat un- 
usual and a magnified view of this region of the phase 
diagram is shown in Figure O At the triple point, the 
sc crystal shows a slightly lower density than the liquid, 
the bcc crystal being the denser phase. The sc crystal 
is thermodynamically stable in a very narrow window of 
temperatures above the triple point, and in this region 




P 



FIG. 7: T — p phase diagram of our octahedral six-patch 
particle system (with a =0.3). Labels show the region of 
stability of each phase. The points at which the reentrant 
behaviour occurs are indicated with a cross. 
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FIG. 8: p — T phase diagram of the octahedral six-patch parti- 
cle system (with a =0.3). Labels show the region of stability 
of each phase. The reentrant behaviour occurs at the points 
at which there is a change of sign in the slope of the phase 
boundaries. These points are indicated by a cross. The black 
circle shows the point where inverse melting occurs. 



the sc-liquid coexistence lines adopt a dome-like shape. 
At the top of the dome, the liquid and the sc crystal are 
in equilibrium, both with the same density. The tran- 
sition is first order though, because the enthalpy differ- 
ence between both phases is not zero. For a small range 
of temperatures below this point, there are two values 
of the pressure at which the liquid and sc phases are in 
equilibrium. At the lower coexistence pressure, the sc 
crystal is the more dense phase. However, the situation 
is reversed at the higher coexistence pressure, where the 
liquid is in equilibrium with a lower density sc phase. 
This means that the coexistence curve shows a reentrant 
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FIG. 9: (Colour online) Detailed view of the phase diagram 
in the region of the sc-liquid-bcc triple point. The triple point 
is shown with a dashed line. Labels indicate the region of sta- 
bility of each phase. The point at which reentrant behaviour 
occurs is indicated by a cross. 
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FIG. 10: Equation of state for the liquid phase (circles) and 
the sc crystal (asterisks) at T* = 0.2. 



behaviour in the vicinity of the triple point, i.e. there is 
a change in the sign of the slope of dp/dT (Fig. [SJ at 
the value of the pressure that coincides with the top of 
the dome in the T — p diagram (Fig. [9]). This kind of 
reentrant behaviour, first speculated by Tammann to oc- 
cur for water (see Ref^), has been found by computer 
simulations both for primitive models of water™ and for 
realistic models of water.— The origin of this behaviour 
can be seen in Figure ITTJl Since the compressibility of the 
solid is very small, the fluid can be less dense and more 
dense than the solid at different points along an isotherm. 
This leads to the existence of two coexistence pressures 
for a given temperature (one low and one high), which 
results in a reentrant behaviour. The low compressibility 



Phase 1 Phase 2 Phase 3 T* p* pi p* 2 p% 
sc liquid bcc 0.229 0.758 0.709 0.728 1.110 

liquid fcc-d (PC) bcc 0.336 3.42 0.956 1.045 1.125 



TABLE VII: Thermodynamic states for the two triple points 
in the phase diagram. 

of the solid is related to the strong directionality of the 
patchy bonds in our model (or hydrogen bonds in the 
case of the water models). 

At a higher temperature, there is another triple point, 
where the fluid, the fcc-d plastic and the bcc crystal are 
at equilibrium. Curiously, the fcc-d plastic phase can be 
less dense than the bcc crystal. As the plastic phase is 
favoured by its high entropy not energy (Table |TT|, the 
density can change considerably without any accompa- 
nying large changes in the energy (as long as the repul- 
sive cores of the particles do not overlap). By contrast, 
the bcc phase is stabilized by its low energy, and devia- 
tions of the second nearest neighbour distance from the 
minimum in the LJ potential will result in a significant 
energetic penalty. As a consequence of the wider range 
of densities possible for the disordered fcc-d phase, the 
coexistence curves again show a reentrant behaviour and 
adopt a dome-like shape in the p—T phase diagram. 

It is also worth mentioning that the bcc-fcc-o coex- 
istence line seems to show inverse melting at very low 
temperatures (T* w 0.08, see black circle in Fig. [H]). 
This inverse melting is different from the reentrant melt- 
ing mentioned above in the sense that inverse melting is 
not caused by the fact that both phases have the same 
density (volume), but because they have the same en- 
thalpy (see Eq. 0|) . Tammann also speculated about this 
possibility,''— and this unusual inverse melting behaviour 
has previously been observed for 3 He and A He and for 
poly (4-methylpentene- 1 ) >^i^ 

The results of this work are consistent with previous 
calculations using a somewhat similar six-patch modelfii 
albeit with differences due to the differences in the poten- 
tials used. Chang et al. also found several solid phases, 
including the plastic fcc-d phase at high temperatures, an 
ordered fee solid at moderate pressures and low tempera- 
tures, and a sc solid at low temperatures and pressures^ 

IV. CONCLUSIONS 

We have computed the phase diagram of particles with 
six octahedrally-arranged patches for a case where the 
patches are relatively narrow. Even for this relatively 
simple potential, a complex phase diagram results with 
competition between multiple solid phases differing both 
in their density and whether they are orientationally or- 
dered, and with unusual features such as reentrance and 
inverse melting. 

Consistent with previous results ; 7 ' 10 ^ 1 there is no sta- 
ble gas-liquid phase separation in our model. The lowest- 
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energy structure is a simple cubic crystal, but this struc- 
ture is only stable at relatively low pressure. As the pres- 
sure increases, first a bcc and then an fee phase become 
most stable, where both phases are orientationally or- 
dered. However, on increasing the temperature, entropic 
effects become more important and a plastic fee crystal 
becomes most stable. 

Even though we have considered very simple 
anisotropic models, they are able to predict the formation 
of low density crystals, in analogy with the preference of 
the proteins to form open crystals. In that sense, even 
though our model is still far from real proteins, our work 
can potentially provide some insights into the complex 
fluid-solid equilibrium of proteins. In future work, it will 
be interesting to explore how the geometry and number 
of the patches will influence the phase diagram, and in 
particular the structure of the stable crystalline phases. 
Particularly relevant to proteins may be the exploration 
of random, rather than just ordered, arrangements of the 
patches. 

Although our calculations have determined the region 
of thermodynamic stability for each phase, it does not 
necessarily follow that these phases will be easily accessi- 
ble from within these regions. Indeed there is increasing 
evidence that the dynamics of crystal nucleation can de- 
pend sensitively on the nature of the crystalline phase 
and also of the liquid ) 66 ' 67 ! 68 ' 69 as has been seen in pre- 
liminary calculations for the current models Therefore, 
in future work we are planning to study the nucleation 



dynamics for the current model, and in particular to ex- 
plore how this dynamics depends on the geometry of the 
patches. Such information might provide important in- 
sights that could help colloidal chemists in designing their 
anisotropic particles to crystallize into the desired target 
structure. 
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